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G\ Abstract 
O 

The present obliquity of Mercury is very low (less than 0.1°), which led previous studies to always adopt a nearly zero obliquity 
during the planet's past evolution. However, the initial orientation of Mercury's rotation axis is unknown and probably much 
j^Q different than today. As a consequence, we believe that the obliquity could have been significant when the rotation rate of the planet 
^ first encountered spin-orbit resonances. In order to compute the capture probabilities in resonance for any evolutionary scenario, 
<^ we present in full detail the dynamical equations governing the long term evolution of the spin, including the obliquity contribution. 
The secular spin evolution of Mercury results from tidal interactions with the Sun, but also from viscous friction at the core- 
mantle boundary. Here, this effect is also regarded with particular attention. Previous studies show that a liquid core enhances 
drastically the chances of capture in spin-orbit resonances. We confirm these results for null obliquity, but we find that the capture 
1 probability generally decreases as the obliquity increases. We finally show that, when core-mantle friction is combined with 



obliquity evolution, the spin can evolve into some unexpected configurations as the synchronous or the 1/2 spin-orbit resonance. 
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1. Introduction 



The present rota ti on rat e of Mercury was discovered by 
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Pettengill and Dvcd ( 1 19651) . when using the new planetary 
radar at Arecibo Observatory in Puerto Rico. Subsequent 
observations confirmed that contrary to previous expectations 
fechiaparelliL [l890; De francescol Il988h . the rotation of this 
planet was not synchronous with the orbital mean motion, but 



presented a peculiar 3 /2 resonant equilibrium ( Mc Govern et al 



119651: ICoTo mbo. 1965). Within a year of the discovery the sta 
bility of this equilibrium became understood, as the result of 
the solar torque on Mercury's qu adrupolar moment of inertia 
pombined with an ec c entric orbit (Colombo and Shapiro, 1966; 
iGoldreich and P eale. 1966; ICounselman a nd Shapiro, 1970). 
However, the reason why this state initially arose remained un- 
satisfactory for a long time. 

Mercury like all the other planets in our Solar System is sup- 
posed to have had an initially rapid spin, that was slow ed down 
by the continuous action of the intens e solar tides (IDarwin , 
18805 IPealel 1 1974 1 19761: iBurnsl Il976l) . As the spin rate ap- 



remaine d somewhat unsatisfactory. 

Later, ICorreia and Laskarl (120041) have shown that, as the or- 
bital eccentricity of Mercury is chaotically varying, with some 
excursions to high values, the rotation rate of the planet can 
be accelerated again, and the 3/2 resonance could have been 
crossed many times in the past. Performing a statistical study 
of the past evolution of Mercury's orbit, over 1000 cases, it was 
demonstrated that capture into the 3/2 spin-orbit resonant state 
is in fact, and without the need of a specific core-mantle ef- 
fect, the most probable final outcome of the planet's evolution, 
o ccurring about 55.4% of the time. 

Goldreich and Peale ( 1967) had nevertheless pointed out that 
the probability of capture could be greatly enhanced if a planet 
has a molten core. In 1974, the discovery of an intrinsic mag- 
netic field by the Mariner 10 spacecr aft seemed to imply th e 



existence of a conducting molten core (Nes s et all 1 1974. 1975) 



and more recently iMar got et alJ ((2007) confirmed its existence 
by using radar observations. Core-mantle friction is then an 
effect to take into account and we expect an increment in the 
captur e probabilities for the 3/2 r esonance. However, accord- 



proaches the orbital m ean motion, it will cross a s eries of reso- 
nances. In their work, Goldreich and Peale ( 19661) have shown 
that, since the tidal strength depends on the planet's rotation 
rate, it creates an asymmetry in the tidal potential that allows the 
capture into these spin-orbit resonances. They also computed 
the capture probability into these resonances for a single cross- 
ing, and found that for the present eccentricity value of Mercury 
(e = 0.206), and unless one uses an unrealistic tidal model with 
constant torques, the probability of capture into the present 3/2 
spin-orbit resonance is on the low side, at most about 7%, which 



ing to IGoldreich and Peald (119671) . this also in creases the cap 



ture pr obability in all the previous resonances. IPeale and Boss 
d 19771) indeed remarked that only very specific values of the 
core viscosity allow to avoid the 2/1 resonance and permit cap- 
ture in the 3/2. 

More recently ( Correia and Laskai . 20091) . it was shown that, 
as the chaotic evolution of Mercury's orbit can also drive its 
eccentricity to very low values during the planet's history, any 
previous capture can be destabilized whenever the eccentricity 
becomes lower than a critical value, except for the 1/1 reso- 
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nance. Including the core-mantle friction effect combined with 
the chaotic evolution of the eccentricity, it was found that the 
spin ends 99.8% of the time captured in a spin-orbit resonance, 
mainly distributed by the following three configurations: 5/2 
(22%), 2/1 (32%) and 3/2 (26%). Although in this case the 
present 3/2 spin-orbit resonance is not the most probable out- 
come, it was also shown that the probability of ending up in 
this resonance can be increased up to 55% or 73%, if the eccen- 
tricity of Mercury in the past has descended below the critical 
values 0.025 or 0.005, respectively. 

The present paper is the conti nuation of the previous works 
( Correia and Laskai , 2004 , 20091) . Here, we revisit in full detail 
the theory of dissipative effects (tides and core-mantle friction) 
and mechanisms for capture in resonance, in order to compare 
with the results from previous stu dies, in particular those from 
Goldreich and Peale (1966, 1967). We are part cularly inter- 
ested in considering the effect of a non-zero obliquity, since all 
previous studies assumed that the obliquity was close to zero, 
because this corresponds to the final evolution from dissipative 
effects. However, at the time of the first resonance crossing, the 
obliquity may still be significant, which can lead to important 
modifications in capture probabilities. In a forthcoming paper 
we will present the full dynamics of the spin with planetary per- 
turbations, which are not included in the present work. 

In the next section, we give the averaged conservative equa- 
tions in a suitable form for simulations of the long-term vari- 
ations of Mercury's spin, including the resonant terms and the 
precession motion. In Section 3 we define a model for tak- 
ing into account the tidal effects including the obliquity con- 
tribution. Section 4 is devoted to the analysis of the core- 
mantle friction. This effect can be divided in two parts, one 
resulting from the libration around the resonance (included by 
iGoldreich and Peale] dl967lV ). and a non-resonant term which 
depends on the obliquity. While the first term tends to increase 
the chances of capture, the other has the opposite effect. In Sec- 
tion 5 we revisit the theory of spin-orbit resonances and Section 
6 is devoted to dynamical equation analysis and its implications. 
In Section 7 we perform some numerical integrations, tracking 
the spin evolution from its origin to the present and computing 
the chances of capture in different resonances. 

2. Conservative motion 

We will first omit the dissipative effects, and describe the spin 
motion of the planet in a conservative framework, including the 
obliquity contributions. The motion equations will be easily 
obtained from the Hamiltonian function of the total rotational 
energy of the planet. 

Mercury is considered here as an homogeneous rigid body 
with mass m and moments of inertia A < B < C, supported 
by the reference frame (i, j, k), fixed with respect to the planet's 
figure. We do the gyroscopic approximation, i.e., we merge the 
axis of principal inertia and the axis of rotation, since for a long- 
term study we are not interested in nutations. Let L be the total 
rotational angular momentum and (I, J, K) a reference frame 
linked to the orbital plane (where K is the normal to this plane). 
The angle between k and K is the obliquity, e, and thus, cos e = 



k ■ K. The Hamiltonian of the motion can be written using 
canonical And oyer's action var i ables (L, X) and t heir conjugate 
angles (6,-iff) dAndoverl [l92l iKinoshitai 1 19771) . L = L k = 
Ceo is the projection of the angular momentum on the C axis, 
with rotation rate a> — 8 - iff cos s and X = L ■ K its projection 
on the normal to the ecliptic; 8 is the hour angle between the 
equinox of date and a fixed point of the equator, and iff is the 
general precession angle (FigJT}. 
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Figure 1: Andoyer's canonical variables. L is the projection of the total rota- 
tional angular momentum vector L on the principal axis of inertia k and X its 
projection on the normal to the ecliptic K. The angle between the equinox of 
date y and a fixed point of the equator A is the hour angle 0, and i/i = y N + N 
yo is the general precession angle. The direction of yo is on a fixed plane E c q, 
while y is on the mean orbital (or ecliptic) E ct of date 



2.1. Gravitational potential 

The gravitational potential *V generated b y the planet at a 



generic poin t of the space r is given by (e.g. iTisserandl 11891 
SmartLll953h : 
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(1) 



where G is the gravitational constant, r — r/r and P2W = 
(3x 2 - l)/2 are the Legendre polynomials of degree two. The 
potential energy 1A when orbiting a central star of mass m Q is 
then: 

<U = m Q <V{r) . (2) 

For a planet evolving in a non-perturbed keplerian orbit, we 
write: 

r = cos(?zr + v)I + sin(tzr + v)J , (3) 

where m is the longitude of the perihelion and v the true 
anomaly. Thus, transforming the body equatorial frame (i, j, k) 
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in the ecliptic one (I, J, K), we obtain (e.g. Correial 2006 ): 
( f ■ j = — cos w sin + sin w cos cos e , 
r ■ k = - sin w sin e , 



(4) 



where w = m+if/+v is the true longitude of date. The expression 
for the potential energy (fJJ becomes: 
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where 



F(0, w, s) = 2 cos(20 - 2w) cos — 



+2 cos(20 + 2w) sin 4 ( |J + cos(26>) sin 2 e 



and 



C 



i(A + B) 



k/R 5 



u> + 5E d 



(6) 



(7) 



C 3GC 
R is the planet's radius and kf the fluid Love number (pertaining 
to a perfectly fluid body with the same mass distribution as the 
actual planet). E d is the dynamical ellipticity, the first part of 
this expressio n corresponding to the flattening in hydrostatic 
equilibrium dLambecM fl98ol) . and 5Ed to the departure from 
this equilibrium. 

2.2. Averaged potential 

Since we are only interested in the study of the long-term 
motion, we will average the potential energy H over the rota- 
tion angle and the mean anomaly M, after expanding the true 
anomaly v in series of the eccentricity e and mean anomaly. 
However, when the rotation frequency a> 6 and the mean mo- 
tion n — M are close to resonance {lo pn, for a semi-integeiQ 
value p), we must retain the terms with argument 2(0 - pM) in 
the expansions 



cos(20) 1 



= — ■ V G(p, e) cos(26» - 2pM) (8) 
a i t—L 

p=—oo 



and 



cos(26> - 2v) 1 



= — Y H(p,e)cos(26-2pM) , (9) 

where a is the semi-major axis of the planet's orbit and the func- 
tions G(p, e) and H(p, e) are power series in e (Tab. [TJ. The 
averaged potential 14 becomes: 
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cjx 2 B 



(1 - x 2 ) G(p, e) cos 2(0 - pM) 



We have retained the use of semi-integers for better comparison with pre- 
vious results 



+ H(p, e) cos 2(0 -pM-< 
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■ H(-p, e) cos 2(0 - pM + 4>) 



where x-Xjh- cos s, (p — us + t//, 
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E d 
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is the 'precession constant' and 

3Gm B —A 
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We can rewrite expression ( TTOb simplified as: 

ni ojx 2 b x 
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(10) 



(11) 



(12) 



(13) 



where the amplitude B x and the phase angle <p x are functions 
depending on both x and ifr, whose expressions are given in ap- 
pendix [A] 

2.3. Equations of motion 

The Andoyer variables (L, 0) and (X, —if/) are canonically 
conjugated and thus 

dL__daj_ <ix__MA_ d4_ _ _m_ 

~dt~~~d0' ~dt ~ ~b\f7 ' 1t~~~dX' 

Despite their practical use, Andoyer's variables do not give a 
clear view of the obliquity variations. Since cose = XfL they 
can be obtained as: 



. ds _ 1 IXdL dX\ _ 1 
dt L\L dt dt) L 

Then, from equation (TTOb we get: 



X 30 + dif/ 



and 



dh 

— = -CB X sin 2(0 - pM - <p x ) , 
dt 

ds 

— = -a r sin s cos 2(0 - pM - <p r ) 
dt 

dif/ 

— = a x + a r sin 2(0 — pM - (f> r ) , 
dt 



(15) 

(16) 
(17) 

(18) 



where a r and (f> r are functions depending on both x and if/, 
whose expressions are given in appendix |A] For non-resonant 
motion the previous equations simplify as: 

dL ds dif/ 

— = — =0 and — = a cose. (19) 
dt dt dt 

The planet spin motion reduces to the precession of the spin 
vector about the normal to the orbital plane with rate a cos e. In 
general we have a » a r and the precession rate including the 
resonant motion ("EqJTSb is nearly the same as the non-resonant 
case (Eq JT9b. 
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Table 1: Coefficients of G(p,e) and H(p,e) to e 4 . The exact expression of these coefficients is given by G(p,e) = ^ £ (j) exp(i2pM)dM and H(p,e) = 
1 ( f ) 3 exp(i 2v) exp(i 2pM) rfM. 



3. Tidal effects 

Tidal effects arise from differential and inelastic deforma- 
tions of the planet due to the gravitational effect of a perturbing 
body. Their contributions to the spin variations are based on 
a very general formulation of the tidal potential, initiated by 
George H. Darwin (1880). The attraction of a body with mass 
m Q at a distance r from the center of mass of the planet can be 
expressed as the gradient of a scalar potential 'V', which is a 
sum of Legendre polynomials: 



(20) 



1=2 1=2 v 7 



where r' is the radial distance from the planet's center, and S 
the angle between r and r' . The distortion of the planet by this 
potential gives rise to a tidal potential, 



2>f 



(21) 



1=2 



where e Vf = fc/'VJ at the planet's surface and k/ is the Love num- 
ber for potential. Since the tidal potential *V? is an Zth degree 
harmonic, exterior to the planet it must be proportional to r 
(solution of a Dirichlet problem). Furthermore, as upon the sur- 
face r' — R <sc r, we can retain in the expansion only its first 
term, 1-2: 

T « = _ t! £^(£) 3 (£) J p! ,c„ sS >. (2 2) 



In general, imperfect elasticit y will caus e the phase angle of 
<V g to lag behind that of 'V' dKaulal 1 1964 by an angle 5{cr) 
such that: 

o-Af(o-) 

5(cr) = — - — , (23) 

Af(cr) being the time lag associated to the tidal frequency cr (a 
linear combination of the inertial rotation rate a> and the mean 
orbital motion n). 

3.1. Equations of motion 

Expressing the tidal potential given by expression d22l) in 
terms of Andoyer angles (6, iff), we then easily obtain its contri- 
bution to the spin evolution as: 



dh / d r V g dX ,d r V 8 
— = —m — — ; — = m — — ; 
dt 36 dt b\j/ 



(24) 



where m' is the mass of the interacting body. As we are inter- 
ested here in the study of the secular evolution of the spin, we 
will average (1241) over the periods of mean anomaly, longitude 
of node and perihelion of the perturbing body. When the inter- 
acting body is the same as the perturbing one (m! = m ), we 
obtain: 



dh Gm%R 5 
dt ~ a 6 



2r>5 



o-)® L M, 



0: 



(25) 



Gm%R 



Yjbio^KM , (26) 



cIe 
dt 

whe re the coeff icients &, T (x, e) are polynomials in the eccentric- 
ity ( Kaulal 1964 ). The factors b(cr) are related to the dissipation 



of the mechanical energy of tides in the planet's interior respon- 
sible for the time delay Af(cr) between the position of "maximal 
tide" and the sub-solar point. They are related to the phase lag 
(5(cr) as: 



b(cr) = k 2 sin 2d(cr) = k 2 sin (crAf(cr)) . 



(27) 



Dissipation equations d25i l and ( f26b must be invariant under the 
change (co, s) by (-to, n—e) which imposes that b(cr) = —b(—cr), 
that is, b(cr) is an odd function of cr. Although mathematically 
equivalent, the couples (co, e ) and (-co, n- e) correspon d to two 
different physical situations dCorreia and Laskai , 2001 ). 



3.2. Dissipation models 

The dissipation of the mechanical energy of tides in the 
planet's interior is responsible for the phase lags 6(cr). A com- 
monly used d imensionless measure of tidal d amping is the qual- 
ity factor Q jMunk and MacDonald , I960!), defined as the in- 
verse of the "specific" dissipation and related to the phase lags 
by 

2^= TF = cot2 <^' (28) 
AE 

where E is the total tidal energy stored in the planet, and AE 
the energy dissipated per cycle. We can rewrite ( l27T i as: 



Her) 



k 2 sign(cr) 
V2V) + 1 



sign(cr) 



Q(o-) 



(29) 



The present Q value of the planets in the Solar system can be 
estimated from orbital measurements, but as rheology of the 
planets is badly known, the dependence of b(cr) on the tidal 
frequency cr is subject to various approximations. 

3.2.1. The visco-elastic model 



Darwinl dl908l) assumed that the planet behaves like a 
Maxwell solicH of constant density p, and found: 



b(cr) = kf- — h — — ^-to~ . 



l+(T b of 



(30) 



where r a = v e /fx e and are the time constants for damping of 
the body tides, 



T b = T a (\ + l9fi e R/2Gmp) . 



(31) 



The visco-elastic model is a very realistic approximation of 
the planet's deformation with the tidal frequency. However, 
when replacing expression d30l into the dynamical equations 
d25l > and ( f26b we get an infinite sum of terms. This problem 
can be solved by using simplified versions of the visco-elastic 
model for specific values of the tidal frequency cr. For instance, 
when cr is small, (jb cr) 2 can be neglected in ( 1301 and b(cr) be- 
comes proportional to cr. 



3.2.2. The viscous or linear model 

In this model, it is assumed that the response time delay to the 
perturbation is independent of the tidal frequency, i.e., the po- 
sition of the "maximal ti de" is shifted from the sub-solar point 
by a constant time lag At (IMignardi 1 197911 19801) . As usually we 
have crAf <K 1, this model becomes linear: 



b s (cr) = k 2 sin(crAf) ^ k 2 crAt 



(32) 



Substituting the above formula into expressions (l25T l and (l26l l, 
we simplify the motion equations as ( appendix lETi: 



dL 
dt 



— = -CK 



1 + X CO 

——n(e)-~xN(e) 
2 n 



ds 
dt 



-7- = K- 



sine 



CO 



xQ.(e)— -N(e) 
2n 



(33) 



where 



Q.(e) = 



1 + 3e 2 + 3e 4 /8 



N(e) = 



(\-e 2 f 2 ' 
l + 15e 2 /2+45e 4 /8+5e 6 /16 

%Q\ml\a) 



(34) 
(35) 
(36) 



Q' 1 = nAt and ^ = C/(mR 2 ). The viscous model is a particular 
case of the visco-elastic model and is specially adapted to de- 
scribe the behavior of planets in slow rotating regimes (to ~ n). 



3.2.3. The constant-Q model 

Since for the Earth, Q changes by less than an order of mag- 
nitude between the Chandler wobble period (about 440 days) 
and seismic periods of a few seconds, it is also common to treat 
the specific dissipation as independent of frequency. Thus, 



1/(0-) - sign(cr)£ 2 /e . 



(37) 



For long-term evolutions and slow rotating planets, this model 
is not appropriate as it gives rise to discontinuities for cr = 0. 
However, it can be used for periods of time where the tidal fre- 
quency does not change much, as is the case for fast rotating 
planets. Substituting Eq.(l37l> into expressions d25b and d26l >. we 
can simplify the motion equations as: 



(38) 



dL 




CK , 


dt 


= -h 




de 


K 


sin s 


~di ' 


16 


CO 



where s g = sign(L) = sign(w). 



2 A material is called Maxwell solid when it responds to stresses like a mass- 
less, damped harmonic oscillator. It is characterized by a rigidity (or shear mod- 
ulus) /.f f and by a viscosity v e . A Maxwell solid behaves like an elastic solid 
over short time scales, but flows like a fluid over long periods of time. This 
behavior is also known as elasticoviscosity. 



4. Core-mantle friction effect 

The Mariner 10 flyby of Mercury revealed the presence of 
an intrinsic magnetic field, which is most likely due to motions 



5 



in a conduc t ing flu id inner core (for a review see iNessl Il978t 
Spohn et all 120011) . Subsequent observations made with Earth- 



•e, 



based radar provided strong evidence that the mantle of Mer- 
cury is decoupled fr om a core that is at least partially molten 
(Mar got et all 120071) . If there is slippage between the liquid 
core and the mantle, a second source of dissipation of rotational 
energy results from friction occurring at the core-mantle bound- 
ary. Indeed, because of their different shapes and densities, the 
core and the mantle do not have the same dynamical ellipticit 
and t he two parts tend to precess at different rates (iPoincan 
Il910h . This tendency is more or less counteracted by different 
interactions produced at their interface: the torque N of non- 
radial inertial pressure forces of the mantle over the core pro- 
voked by the non-spherical shape of the interface; the torque 
of the viscous (or turbulent) friction between the core and the 
mantle; the torque of the electromagnetic friction, caused by 
the interaction between electrical currents of the core and the 
bottom of the magnetized mantle. 

4.1. Equations of motion 

We will adopt henceforward a model for the pl anet which is 
an extension of the model from IPoincare ( 1910h of a perfect 
incompressible and homogeneous liquid core with moments of 
inertia A c = B c < C c inside an homogeneous rigid body with 
moments of inertia A m < B m < C m , supported by the same ref- 
erence frame (i,j,k), fixed with respect to the planet's figure 
(FigtD- The combined effects of inertial and frictional cou- 
pling across the ellipsoidal core-mantle boundary are taken into 
account, assuming laminar flow. 

Denoting 6 = ia—co c the differential rotation between the core 
and the mantle, we can write the non-radial inertial pressure 
torque of the mantle on the core in a general f ormulation to first 
order in the core dynamical e llipticity, E c , as dRochesterl[l976k 
Sasao et allll980t|Pais etailll999h : 



N = co c x L c - C C E C wkxtf. 



(39) 



where L c - I c ■ a> c is the core angular momentum with I c - 
diag(A c , A c , C c ) its tensor of inertia. 

The two types of friction torques (viscous and electromag- 
netic) depend on the differential rotation between the core and 
the mantle and can be expressed by a single effective friction 
torque, T. As a general expression for this torque we adop t 
dRochesteii[l976tlsasao et1ulll980kliylathews and Guol 12005) 



T = C c (k + k' kx) S , 



(40) 



where k and k' are effective coupling parameters. They re- 
sult either from viscous and electromagnetic stresses at the 
core-mantle interface and can be written as a sum of these 
two effects: k = k v \ s + K em and k' - k' ■ + ic' em . Recent es- 
timations of these coeffi cient s can be found in the work s of 
Mathews and Guol d2005l ) and beleplace and Cardinl d2006l) . In 
the simplified case of no magnetic field, the coupling parame- 
ters are only given by th e viscous f r iction c ontributions, which 



can b e simplified as ( Noi r et al 
2005): 



2003; Mathews and Guo 



2.62 y/v\t7\/R c and < is = 0.259 y/v\t7\/R c , (41) 



where R c is the core radius and v the kinematic viscosity, which 
is poorly known. Even in the case of the Earth, the uncertainty 
in v c overs about 13 orders of magnitude dLumb and Aldr idge. 
1 19911) . the best estimate so far being v ^ 10 6 m 2 s 1 flGansl 
1 197a IPoirierl 1 1 988k Iwijs etailll998l) . 

Since the derivative of the angular momentum is given by 
the sum of external torques, the contribution of the core-mantle 
friction is the solution of the system: 



' dL„ 
dt 



= p + t-n-t , 



dL c 

— £ = N + T , 
dt 



(42) 



where P = rx VI/ is the precession torque and T — r x m'W g 
the tidal torque. L m denotes the mantle's angular momentum, 

L m = C m co = C m cok (43) 

and the total angular momentum variations are given by: 

dL dL m dL r 

— = + =P + T . (44) 

dt dt dt 

L precesses around K, a normal vector to the orbital plane, with 
angular velocity £1 given by: 



n = -ij/K + sp , 



(45) 



where e accounts for the secular effects resulting from tides and 
core-mantle friction and p = Kxk / sin s is the unit vector along 
the direction of the averaged precession. Thus, 



dL 
~dt ~ 

Projecting it over 



Q x L = il x (Ccok -Ic -8) 



q = kxp 



K - kcos s 



sine 



(46) 



(47) 



we get 

(a XL) ■ (kxp) = (£2 • k)(L ■ p) - (£1 ■ p)(L ■ k) 

= -ifi cos e{—p I c S) — e (Cw - k ■ I c ■ 6) 

= ij/ cos sA c 6 p - s(Cu> - C c 6k) , (48) 

with (5k = 6 ■ k and 5 P — 6 ■ p. Then, assuming ^ <K w we have 
from expression d44l i 



P q T q ij/ cos eA c 5 p 
Cto Cto Ca> 



(49) 



where —PJCto and -Tq/Cto are respectively given by expres- 
sions (TTTb and d26l ). T he value of 6„ in the case of uniform p re- 
cession is ( Rochester, [l976l; IPais et all 1 1 999k ICorreial [2006 ) : 



Ktoijj sin s 



p ~ .,2 



K + (k + E c L0) 



(50) 



In the case of a non-uniform precession, as it can be the case of 
Mercury close to a spin-orbit resonance. ICorreial (120061) shows 
that the asymmetric terms in (B - A) are periodic and average to 
zero, and we can still use the previous expression for 5 p . Using 
c c = C c /C, A c =i C c and ip acose, we finally have for the 
obliquity variations with core-mantle friction: 



Ceo Ceo 



Kf cos 3 e sin e . 



with 



c c Ka 



k 2 + (k + E c Cof 



(51) 



(52) 



which is always a positive quantity. 

The rotation rate variations of the mantle and the core can be 
obtained by projecting both equations d42b onto k: 



inside the integral of previous expression only i\ is not secular, 
we have 



1 

«5k(«) = - 

K 



Pk C m . 

ilfsmsdr, 

C C 



+ P(t), (59) 



with 



-I 



Pk 



P(t) = &-*■»' | — e"'"'dt. 

Cm 



(60) 



4.2.1. Strong coupling 

We consider a strong coupling between the core and the man- 
tle whenever k 2 » B, where B is given by expression ( fT2b . In 
this case, we can simplify expression d60b by performing an 
integration by parts: 



d dhs afk dLi 
—(Li • k) = k + L, — = 

dt dt dt dt 



k + Li ■ (« x k) . (53) 



For the mantle, L, is given by expression d43l and then 



dco 

C m — = Pk + T k - C c k 6 k , 
dt 



(54) 



where Pk and 7\ are respectively given by expressions ([Tol l and 
( 1251 ). For the core L, = L c — I c ■ co c and we have 



1 ~ k C k 2 C 2 dt 



(61) 



where we neglected terms higher than B/k 2 . Substituting the 
above equation (loTt into expression (l54l we find for the rotation 
rate variations: 



do P k T k c c 8 P k 

— = 1 o)Kf cos £ sin" £ H . (62) 

dt C C 1 k dt C 



(55) 



daf c 

C c = C c k 5 k + L c ■ {-ib sin sp - e q) 

dt 

— C c k 6 k + A c i}f sin £ 5 P + A c s 6 q , 



where a>^ - a> c ■ k = to - 6 k and 6 P is given by expression i 
We can also neglect the term in A c e 5 q because according to 
expression (l49l its average is a second order term in 5. Thus, 
using A c C c we have 



4.2.2. Weak coupling 

For weak coupling we assume k 2 B, we will thus neglect 
second order terms in k 2 IB. Performing again an integration by 
parts in expression d60l l. but changing the roles of P k /C m and 
s" mt we have: 



Pit) 



dt . 



(63) 



da>: 



— = k 6 k cos 2 £ sin 2 £ . 

dt c c 



(56) 



which gives for the rotation rate variations when substituted into 
expression i 



4.2. Differential rotation 

Combining equations (l54l and (IB3I ) we find a differential 
equation for (5k, 



dco Pk T k 77 
— = — + coKf cos £ sin £ - c, 

dt (^J ft) C 



r 



dt . (64) 



d6 k P k Pk A c . . 

= -K m dk H + iff sin £ dp , (57) 

dt C m Cm Cr 



5. Spin-orbit resonances 



where k„, = KC/C m . Its solution allows us completely to deter- 
mine the spin of the mantle (Eql54b without needing the core 
variations (EqBBIi: 



6t(t) 



(Pk Pk Ac. 
1 tlfsmeSn 

c c c 



'dt . 



(58) 



The secular variations to the spin can be seen as constant for 
short-periods of time. Thus, since among all the contributions 



The resonant equilibrium was first observed in the Moon , 
that i s locked in a 1/1 spin-orbit resonance (e.g. iGoldreichl 
19661) . That other spin-orbit resonances were possible was not 



realized bef ore the discovery of the 3 /2 spin-orbit resonance 



of Mercury (Pettenaill and Dvce, 1965), gi ving rise to several 



3/2 
5). 



detailed studies dColombolfl965l:IColombo and Shap irolll966l : 
Goldreich and PealeT 19661: ICounselman and Shapirol 1970h . 
Such non-synchronous spin-orbit resonances require a large or- 
bital eccentricity, but also, as we will see, low obliquity. 



7 



5.1. Effect on the rotation rate 

Neglecting by now the dissipative effects resulting from tides 
and core-mantle friction, when combining expressions ( TToT l and 
d43t . we can write near a generic spin-orbit resonance {co pn): 



dco 

— = -R m sin 2(0 - pM - , 
at 



v) 



(65) 



with /3 m - /3 x /c m , where /3 X and <f> x are given in appendixlAl and 
c m - C m /C. Let us denotey = 6—pM-<p x . Since co = 8-tfrcos s 
we have y = co — pn — if/ x , with \j/ x = <p x + ip cos s. Because we 
assume \p a cos s (EqfT9ll. for small variations of co and e, we 
can consider /3 X , <j> x and ifr x as constants. Thus, we have y - co 
and expression d65l > can then be rewritten as 



y + f3 m sin 2y - , 



(66) 



which is the same as the equation of a free pendulum (Fig|2j. 
The first integral of this equation is given by 



h-y- f3 m cos 2y . 



(67) 



where h is a constant of the motion related to the energy. The 
separatrix equation is given by h - /3 m , where h > /3 m gives 
the trajectories in the circulation zone (outside the resonance) 
and h < /?,„ the trajectories in the libration zone (captured in 
resonance). The maximal and minimal libration width, w_ < 
co < co+, are obtained from the separatrix equation (h =/?„,): 

co± - pn± Aco with Aw = ~<j2/3 m , (68) 

Since < j3/c m (Eq ll241 i. from expression (fT2l) we have: 

Aco 



< l 3 LA*0.02. 

ft V ^ m 



(69) 



using (B-A)IC m a 1.2 X 1(T 4 dAnderson et al.L ll987). 



For the tidal torque we can use the viscous model approx- 
imation (section 13.2.21 ). since significant spin-orbit resonances 
only occur in the slow rotation regime (co ~ n). This torque 
is then given by expression (l33l >. which can be rewritten using 
co - pn + y and E(e) = N{e)/Q(e) as 



T k \+x l 
— = -K 

C 2 



( 2xE{e) 


7 


[ l+x 2 ] 


+ - 

n 



(72) 



According to expression ( 1691 ) we have y/n <s 1. Thus, the 
non-resonant core-mantle friction torque can also be made lin- 
ear: 



coKf(co) =i pnKf(pn) 



J 

1+%- 



where 



Q.f = n 



\n{coK f ) 



tt)=pn 



£ 
P 



(73) 



(74) 



and q is a semi-integer like p. Indeed, since a 2 cc hT 1 (EqQj} 
and k cc co l l 2 (EqHTt. from expression (l52l we have for weak 
friction (k <k E c co) that 



coKr oc — , 

J (xfl 



(75) 



with q — 5/2 and for strong friction (k » E c co) the same previ 
ous expression, but q - 3/2. This result also wo rks for turbu 
lent friction, for wh ich k cc 1 /co and thus q = 4 (lYoderi 1 1995 



Correia et al., 2003) 



Finally, for the resonant core-mantle friction contribution we 
will split our analysis for strong and weak coupling as in sec- 
tion |4j2] In the first situation (k 2 » p) this calculus is trivial 
when using expression d65l l for the precession torque. Indeed, 
from expression doTT i we have 



5.2. Dissipative torques 

Spin-orbit resonant configurations result from an evolution- 
ary process. It is believed that the terrestrial planets' rotation 
was faster at the time of their formation, but due to dissipative 
torques it has decreased (section loTTT l and may have been cap- 
tured inside a resonance when crossing it. 

The secular variations of the rotation rate are easily computed 
from the mantle's angular momentum variations (Eqs J54l59l as: 



where 



dco Pk — 
— = — +D . 

dt C m 



D = - coK f {\ - x 2 )x 2 - c c K m P(t) 



(70) 



(71) 



denotes the mean dissipative torque (x = cos s), which is com- 
posed of three terms: the first arising from tidal effects, the 
second from non-resonant core-mantle friction and the last one 
is a dissipative term resulting from the presence of spin-orbit 
resonances and core-mantle friction together. 



CcKmP(f) 



-c c — + 

C m K 0t 



d (1\ 

c 



c c p ln sin 2y y cos 2y 

K 



(76) 



The term c tV 6,„ sin 2y can be summed with the initial precession 
torque (Eqj65l providing a single term with amplitude /3 X = 
p m — c c Pm- Thus, in the case of strong coupling the spin behaves 
as if there was almost no differentiated internal structure, with 
only a small perturbation resulting from the term in cos 2y. 

In the case of weak coupling (k 2 <k /?), if we neglect the 
secular effects (which is possible for short term variations), we 
can write using equation (l66l l 

P(t) = J ^Ldt^ Jydt = y-y , (77) 

where yo is an integration constant. This approximation is valid 
as long as the time interval Af for which we perform the a bove 
integration verifies Af <K 1 /tc m dCorreia and Laskai , 2009). 
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5.3. Capture probabilities 

The total variation of the rotation rate when dissipative 
torques are included is then 



y + fi m sin 2y = D(y) . 



(78) 



The spin-orbit term /3, n sin 2y is commonly known as the 
restoration torque as it will counterbalance the dissipative 
torque D preventing the planet from escaping the resonant con- 
figuration; 

Goldreich and Peale ( 19661) computed a first estimation of 
the capture probability P cap , and subsequent more detailed stud- 
ies prove d their expressi on to be essentially correct (for a re 
view, see 



Henrard, 



1991 . We consider here the planet's orbit 



as a fixed ellipse at the moment it approaches the resonance, 
since the perturbations of the orbital parameters during this 
short period of time do not change the behavior of the planet 
dGoldreich and PealeL[l966l) . 

Differentiating equation (l67l > and replacing it in expression 
d78l ). we obtain 

dh ~-^r-\dy 



dt ^ dt 



(79) 



The "energy" variation of the planet after a cycle around the 
resonance is then given by the function: 



Ah(yi) 



L dt dt + J h 



dh 
dt 



■ dt 



(80) 



"T3 



D(y)dy + 2 D(y)dy, 
ri -Jyi 

where y\ is the y value when the planet crosses the separatrix 
between the circulation and the libration zones, i.e., the y value 
for h — yS m . y2 and 73 are the first two following y values cor- 
responding to y = 0. fi, ?2 and ?3 are the instants of time where 
the previous events occurred, respectively. 

In order to be captured, after a cycle inside the resonance, 
the planet must remain within the libration zone. When de- 
spinning from faster rotation rates, this means that the total "en- 
ergy" of the planet after a cycle, h{yi), must be smaller than the 
separatrix "energy", h—fi m (Fig (2j . Thus, 

h(y 3 ) = h( 7l ) + M(n) =/3 m + M(n) <p m (81) 

and the capture condition becomes: 

M( ri ) < . (82) 

Assuming that the y\ values comprised between -jt/2 and 
n/2 are distributed uniformly, capture inside the resonance will 
occur for "energy" variations within 



{£ : Ah(-n/2) < Ah(yi) < 0} 
from a total of possibilities: 

{fi r : M(-tt/2) < M(n) < Ah(n/2)} 

that is, 

dyi 



M(-|) 



cap 



dyi 



M(-|)-M(f) 



(83) 



(84) 



(85) 



(a) 




1 Y 






















Figure 2: Capture in the y = 6 — pM - <p x resonance in the phase space (y,y). 
In the first situation (a) the planet is captured, while in the second one (b) it 
manages to cross the resonance without being trapped. 



Usually tidal torques are very weak, i.e., \D(y)\ <K \/3,„\, which 
gives y3 ^ -7t/2, y^ — n/2 and h ^ j3 m . From expression d67] ) 
we then write 

y = sign(y) y[2fr n cos y . (86) 



Replacing it in expression 
probability: 



we compute for the capture 



P + 

cap 



D(y = - cos y)<iy 



(87) 



2/3 m cos yjdy 



Following an identical reasoning, we can obtain the expres- 
sion for the capture probability when the spin is increasing from 
lower rotation rates: 



cap 



i + — 



+ V^cosy) dy 



(88) 



2j3 m cos yj dy 



Let us notice that for even torques in y, the capture never 
occurs, while for odd torques it is unavoidable. The capture 
probability must lie between and 1, but often the results given 
by expressions d87b and ( f88l> are outside this interval. In those 
cases, if P cap < then P cap = 0, and if P, 
For a general dissipation torque in the form 



ca p > 1 then P cap = 1. 



D(y) 



-K 



y 

V + (jui + p.7 cos 2y) — 
n 



(89) 



where K, V, p\ and p2 are constants, we compute from expres- 
sions (|87} and (l88l i: 



P ± -2 

1 cap 



n n V 
2 Aa> p 



with M = M1 + ^. (90) 



6. Dynamical evolution 

In this section we analyze the dynamical equations obtained 
in previous sections. The main goal is to describe both evolu- 
tion and final stages for the spin under the effect of dissipative 
effects. 

6.1. Rotation rate evolution 

The secular variations of the rotation rate are given by ex- 
pression ( TTTT ). For a fast rotating planet, the spin is far from any 
spin-orbit resonance and we can retain only the secular dissi- 
pative terms, because i\ = P(t) = 0. In this regime we can 
use a constant- Q model as good approximation for tidal effects 
(Eql38ll. Since all secular terms involved are negative (with 
to > 0), they can only decrease the rotation rate for any value 
of the obliquity and eccentricity. This is valid until the slow ro- 
tation regime is attained (to ~ n), where tidal effects can coun- 
terbalance the braking effect from core-mantle friction. Once 
in this regime, two different behaviors are possible: the rotation 
rate stabilizes around an equilibrium point given by the solu- 
tion of 7\ = (section or the rotation rate is captured in a 
spin-orbit resonance (section[5]l. 

6.2. Obliquity evolution 

6.2.1. Effect of core-mantle friction 

According to expression (|5T1 >. the secular effect of core- 
mantle friction on the obliquity is given by: 



de 

— = —Kr cos e sin e . 

dt 1 



(91) 



Since Kj > (Eql52t, for any rotation rate the core-mantle 
friction always brings the equatorial plane of the planet to the 
same plane as the orbit. We have that e vanishes for e = 0° and 
e = 180° (which correspond to stable equilibrium positions) 
and for e = 90° (unstable equilibrium). Moreover, since Kj 
is proportional to a 2 (Eql52t and a oc <y _1 (EqlTTV the mag- 
nitude of both to and e will grow as the planet slows down. 
Thus, for fast rotating planets the core-mantle friction effect 
can be neglected, but as the planet arrives in the slow rotating 
regime (co ~ n), this effect grows s o much that it may control 



the entire evolution of the obliquity dCorreia and Laskan 12003 



Correia et al., 2003). The decrease of the rotation rate and the 



obliquity variations are intimately coupled. Indeed, combining 
the secular core-mantle friction contributions from expressions 
( 1731 and d9"Tl i, we have that to cose = we sine (the spin nor- 
mal component is conserved). As a consequence, for an initial 
rotation rate a>j and obliquity e,- + 90°: 



LOi 



COS Si 

cose 



(92) 



Since the obliquity evolves towards 0° or 180°, i.e., | cos e| — * 1, 
the equilibrium rotation rate is attained for to e = to/ | cos e,|. 



6.2.2. Effect of tides 

For a fast rotating planet, using again a constant- Q model, 
the tidal effects on the obliquity are given by the second equa- 
tion in system ( f38b . This equation has a single stable point for 
e = 67.1 1° (or x = 0.388953). Thus, since the core-mantle fric- 
tion effect can be neglected in the fast rotating regime (Eql9TT>. 
whatever is the initial obliquity, it will evolve by tidal effect 
toward this balance point. 

Once the planet arrives in the slow rotation regime (to ~ «), 
the constant- Q model is no longer suitable and expression ( f38T > 
no longer valid. Using the viscous model instead (Eql33l, we 
find that the obliquity still has only one stable point, obtained 
as the solution of ds/dt = 0: 



e = arccos (2nE(e)/to) if to > 2nE(e) . 







if 2nE(e) > to > 



(93) 



where E(e) = N(e)/£l(e) > 1 (FigOJ. Contrary to the fast ro- 
tating regime, here the stable point depends on the eccentric- 
ity and on the rotation rate of the planet. The stable value for 
the obliquity decreases as the rotation slows down, stabilizing 
for circular orbits at zero degrees for rotation rates of 2n or 
smaller, i.e., twice the orbital mean motion. It results that tides, 
like core-mantle friction, always finish by bringing the planet's 
equator to the orbital plane. However, while core-mantle fric- 
tion allows retrograde final rotations, tides alone only admit di- 
rect rotations. 

6.3. Equilibrium positions 

True equilibrium positions for the spin result from simulta- 
neous balance points for the rotation rate and obliquity, that is 



to = and e = . 



(94) 



In section 16.11 we saw that core-mantle friction and tidal ef- 
fects both decrease the rotation rate for a fast-rotating planet. It 
is thus impossible to stabilize the spin before the planet arrives 
in the slow rotation regime. Once in this regim e, core-mantle 
frictio n becomes dominant over tidal effects dCorreia et all 
2003b and drives the obliquity into 0° or 180° (EqED- It fol- 
lows then that we must look for stable values of the rotation rate 
(to - 0) when x — + 1 in order to find the spin equilibrium po- 
sitions. For these obliquity values, in absence of the spin-orbit 
resonant term P^, the contribution of the core-mantle friction 
to the rotation rate vanishes (Eql73l. The equilibrium rotation 
rate is then determined solely by the tidal effects (Eql25ll, that 
is, when 

rk = ° ° J]b(o-)&^(±l,e) = 0, (95) 

cr 

or, using the viscous model (Eq 13311. simply 



to e — E(e) n , 



(96) 



which means that the equilib rium rotation rate increases with 
the ec centricity of the planet dGoldreich and PealeL 11966 : 
1 98 If) . This behavior is illustrated in Figj3] 
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Figure 3: Evolution of the equilibrium rotation rate u) e ln = E(e) with the 
eccentricity when e = 0° using the viscous model. As the eccentricity in- 
creases, co e also increases. The tidal effects lead the planet to exact reso- 
nance when the eccentricity is respectively ei/i = 0, £3/2 = 0.284926803 and 
e m = 0.392363112. 



6.4. Capture in resonance with s = 0° 

In section 1531 we saw that s — 0° or e = 180° correspond to 
the only two stable positions for the obliquity. For those obliq- 
uity values, the non-resonant core-mantle friction contribution 
to the rotation rate vanishes (Eql73l and the rotation rate equa- 
tion (170b can be greatly simplified. Indeed, for s — 0° (x — 1), 
we can write 



P - 
y = H(p, e) sin2y + D (y) - c c K m P{t) , 

£ III 



(97) 



in expression (l90i l. we get for the capture probabilities in the p 
resonance: 



p ± ~ 2 

1 cap 



l g ± l( p-E(e) \ 
p p\ 2Aa>/nn j 



(100) 



wit h p — and n = 1 . Thi s expression was first obtained 
by iGoldreich and Peald (11966b . It is straightforward that if 
(B - A)/C m increases (or decreases), the capture in resonance 
also increases (or decreases). Using the prese nt value for Mer- 
cury's inertia moment (B- A)/ C m 1 ,2x 10~ 4 ([Anderson et al 



19871) and e - 0.206, we compute for the 3/2 spin-orbit reso- 
nance P* = 7.73% and P~ = 0. In figure|4]we plotted several 
examples of capture in different resonances as a function of the 
eccentricity when the planet is de-spinning from faster rotation 
rates CP™) and when its spin is increasing from slower values 
CPcap). The most remarkable feature is that the probability 
grows very fast as the eccentricity increases (or decreases for 
P^p), but it suddenly decays to zero. This impor tant behavior 
was also described by IGoldreich and Peak (1966) and it corre- 
sponds to a zone where the equilibrium rotation rate does not 
reach the resonance circulation zone, that is, when 



\a> e — pn\ > Aa> a « , 3H(p,e) 



B-A 



(101) 



Thus, the resonance width, which depends mainly on (B - 
A)/C m , not only contributes for the probability of capture, but 
also determines whenever this capture can occur or not. 



where y — 9- pM - <p and 

D Q {y) = -KQ.(e) 



p - E(e) + ^ 
n 



(98) 



A similar expression is obtained for the case e = 180° (x = 
-1), we just have to replace to by -a> (see section [3~TT ). For a 
circular orbit (<? = 0), we have H(p + 1,0) = (TabO} and 
thus, the only spin-orbit resonance where capture can occur is 
the synchronous resonance (p = 1). Capture in this resonance 
always occurs, because according to expression d96*l > we have 
a> e = n, since £(0) = 1. 

When the orbital eccentricity increases, according to expres- 
sion d96b the equilibrium rotation rate co e increases to a larger 
value than n (FigO. Since the synchronous resonance width 
(Eq.|69])fore = o is 



Aa> a; n 



3(5 -A) 



r 



(99) 



when a> e > n + Acu, capture in this resonance becomes impos- 
sible. The final rotation rate will then be given by u> e unless 
capture in a resonance with p > 1 occurred. 

6.4.1. Absence of core-mantle friction 

In order to simplify, we will first compute the capture prob- 
abilities when there is no contribution from the resonant core- 
mantle friction, i.e., we neglect the term c c K m P(t) in expression 
d97| >. Thus, replacing the expression of the tidal torque (Eq|98l 



6.4.2. Weak coupling 

According to expression d77l ), in the case of a weak coupling 
(k 2 <sc B) we can rewrite expression (|97| > for e = 0° as 



y = H(p, e) sin 2y + D(y) 

Cm 



(102) 



where 



D(y) = D (y) - c c K m (y - y ) 
= -K£l(e) 



p-E(e)-Q— +p- 
n n 



with 



C c K m fl 

q = and p — 1 + q . 



KQie) 



(103) 



(104) 



Just before the rotation rate enters the libration zone, its aver- 
age value is given by < y >= 2Au)/n (Eql86*l). In this regime, 
k <K ~ Aa>, and therefore the core is unable to follow the 
periodic variations in the mantle's rotation rate. Thus, if the 
planet is de-spinning from faster rotation rates, we can adopt 
yo - 2Aco/n:m expression (1103l l. since capture probabilities are 
computed when the rotation rate crosses the separatrix between 
the libration and the circularization zones (iCorreia and Laskar , 
2009). Likewise, if the planet is increasing its spin from slower 
rotation rates, we use yo - -2Aw/7r in expression ( 1 1 031 >- 

In this case, the capture probabilities are also given by ex- 
pression ( 1 100b . with q > 0. Since p = 1 +g > 1, this expression 
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Figure 4: Probability of capture in some spin-orbit resonances for eccentric orbits with e = 0° and (B - A)/C„, = 1.2 X 10~ 4 , when using the viscous model for 
tides and no core-mantle friction. The dashed line corresponds to a planet increasing its spin from slower rotation rates, while the solid line corresponds to a planet 
de-spinning from faster rotation rates. In this last situation, for all resonances but the 1/2, as the eccentricity increases, capture probability in higher-order resonances 
increases. However, it suddenly decays to zero when the equilibrium rotation rate falls outside the resonance width. For the 1/2 resonance, the planet can only be 
captured when increasing its spin from lower values and with low probability (the maximum being about 1.04% for e = 0.165). 



simplifies to P* dp = 2/// [l ± ( ^^^ )]» that is, the capture prob- 
ability is always higher than in absence of the resonant contri- 
bution from core-mantle friction. In figure[5] we plotted several 
examples of capture in different resonances as a function of the 
eccentricity when the planet is de-spinning from faster rotation 
rates (P^ap) an d wnen its spin is increasing from slower values 
(/"cap) using v = 10 _6 m 2 s . When comparing with figure [4] 
(absence of core-mantle friction), we observe t hat the capture 
probabilities largely increase, as predicted by IPeale and B oss 
(11977b . For instance, with the present eccentricity of Mercury 
(e = 0.206) we get P+ p = 100% of capture in the 3/2 reso- 
nance, which contrasts with P* - 7.73% in the absence of 
core-mantle friction. If we used a stronger value of the viscos- 
ity, v = 10 _2 m 2 s _1 , the capture probability in any of the four 
resonances plotted in figures |4]and[5]is always one if the eccen- 
tricity is higher than 0.002 and smaller than 0.5, including for 
the 1/2 resonance. 

6.4.3. Strong coupling 

Following expression (l76l l. in the case of a strong coupling 
(k 2 » /?) we can rewrite expression (l97T i for s — 0° as 



where 



y = -f3H(p, e) sin 2y + D(y) 



(105) 



D(y) = Do(y) H(p,e)ycos2y 

K 



and 



-KQ.(e) 



E(e) + (1 +jU 2 cos2y)- 
n 



2/3 H(p, e)n 



,(106) 



(107) 



K K£l(e) 

Thus, according to expression (|89l l, the capture probabilities 
are still given by expression (1100b with q — 0, C,„ = C, and 
fi = 1 + jt/2/3. Since p > 1, the consequences for the capture 
probability in resonance are the same as in the weak coupling 
situation (section l6.4.21 i. 

6.5. Capture in resonance with e + 

The final possible evolutions for the obliquity resulting from 
dissipative effects are s = 0° or s = 180° (see section 16731 . 
However, it may happen that when a resonance is crossed, the 
obliquity is still evolving toward one of the final states. It is then 
useful to analyze the consequences of a non-zero obliquity. For 
simplicity, we will first look at the case of a planet without a 
core. The complete effect with core-mantle friction is discussed 
in section [6". 5. 41 

6.5.1. Absence of core-mantle friction 

In absence of the core-mantle friction effect, for a non-zero 
obliquity value, the equilibrium rotation rate u> e is obtained 
from expression ( 1331 ) setting dL/dt = 0: 



2x 



1 + x 



- E(e) n < E(e) n 



(108) 
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Figure 5: Probability of capture in some spin-orbit resonances for eccentric orbits with £ = 0° and (B - A)/C m = 1 .2 X 1CT 4 , when using the viscous model for tides 
and v = 1(T 6 m 2 s~'. The dashed line corresponds to a planet increasing its spin from slower rotation rates, while the solid line corresponds to a planet de-spinning 
from faster rotation rates. When comparing with figure [4] (absence of core-mantle friction), we mainly notice that the eccentricity values that provide 100% of 
chances of being captured largely increase. For the 1/2 resonance, the planet can only still be captured when increasing its spin from lower values but now also with 
higher chances (the maximum being about 36% for e = 0.12). 



Thus, the main consequence of an increasing obliquity is to re- 
duce the equilibrium rotation rate (FigO. For zero obliquity, 
the lowest equilibrium was the synchronous motion (aj e = «), 
so this configuration was the last possible evolutionary stage 
for a planet de-spinning from faster spins. Now, if the spin 
of the planet is not captured in the 1/1 resonance, the planet 
may evolve to slower final spin configurations, including the 
1/2 spin-orbit resonance or below. Indeed, for s > 90° the equi- 
librium is fixed at a negative rotation rate, allowing the capture 
in negative resonances as well (p < 0). Notice, however, that in 
this case the equilibrium spin always corresponds to prograde 
rotation, since s > 90° . 

The capture probabilities in all resonances will also vary 
for different obliquities: not only the tidal torque is different 
(Eql3"3l. but also the resonance width will change CEq J 1 22b . 
From equation ( f90T > we compute: 
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Figure 6: Evolution of the equilibrium rotation rate u) e with the obliquity for 
fixed eccentricities and using the viscous tidal model. The equilibrium rotation 
rate decreases as the obliquity increases. For e = 90° we have u e = for all 
eccentricities. Although we have u) e < for e > 90° , notice that the equilibrium 
spin still corresponds to a prograde rotation state. 



P =2 

1 cap 



p- 



2xE(e) 
1 + x 2 



2Ato 



(109) 



where g = 0, fi = 1, and Au> = 
being given by expression ( 1122l i. 



V2^ = y/2fijc m , with X 
As the obliquity increases, 
the term in x will decrease and the capture probability will then 
also decrease. In figure [7] we plot the evolution of the capture 
probability with the obliquity for the synchronous resonance 
and two different eccentricities when the planet is de-spinning 
from faster rotation rates. For e — (circular orbit) and low 



obliquity the capture probability is 100%. As the obliquity in- 
creases, this probability decreases fast, and it is close to zero 
for obliquities higher than 90°. For eccentric orbits, we observe 
a sort of contrary effect: for some values of the eccentricity, 
the obliquity can be responsible for an augmentation in the cap- 
ture probability. Indeed, for high eccentricity and low obliquity 
the equilibrium rotation rate is always faster than n (Fig(3j) and 
the capture in the synchronous resonance becomes impossible 
(FiglU). As the obliquity increases, the equilibrium rotation rate 
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Figure 7: Probability of capture in the 1/1 spin-orbit resonance when de- 
spinning from faster rotation rates for different obliquities and e = (solid line) 
ore = 0.2 (dashed line). We used the viscous tidal model and (B - A)/C m = 
1.2 x lfr 4 . We see that the capture probability decreases as the obliquity in- 
creases. Therefore, unless £ < 30.66°, the synchronous resonance is no longer 
the last possible stage for the spin evolution. 
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Figure 8: Probability of capture in the 1/2 spin-orbit resonance for eccentric 
orbits when £ = 80°, using the viscous tidal model and (B-A)/C m = 1.2x 10~ 4 . 
The solid line corresponds to a planet de-spinning from faster rotation rates, 
while the dashed line corresponds to a planet increasing its spin from slower 
rotation rates. According to expression 11081 . for £ > 74.5° the equilibrium 
rotation rate oj e is below the 1/2 resonance and for those obliquities P Cd p # 
(Eq[T09t. 



decreases CEq J108b allowing the planet to approach again the 
synchronous rotation rate. For e = 0.2, this occurs when the 
obliquity is around 60° (Figj7j. 



6.5.2. Capture in the 1/2 resonance 

An important consequence of a non-zero obliquity is that 
the spin can reach previously impossible configurations (sec- 
tion l6.51 l. Because the capture in the 1/1 synchronous resonance 
becomes avoidable, the equilibrium rotation rate can evolve to 
lower values than n. In particular, capture in the 1/2 resonance 
can occur. According to expression (1108b . for a planet de- 
spinning from faster rotation rates, unless s > 74.5° the equi- 
librium rotation rate is above the 1/2 resonance and the capture 
probability in this resonance is zero (Eq |1091 >. Thus, there is 
only a chance of capturing the planet if the resonance is crossed 
with high obliquity (Fig IS). 

However, when the initial obliquity of the planet is high, fol- 
lowing expression d92l it can happen that the rotation rate be- 
comes inferior to n/2 when the obliquity is brought to zero de- 
grees by core-mantle friction. Then, the planet rotation rate 
will increase towards the equilibrium given by equation d96l l 
and cross the 1/2 resonance from slower rotations. The capture 
probability is given by P~ ap (Eq llOOt with p - 1/2. Notice also 
that, as for any resonance different from the 1/1, capture in the 
1/2 resonance is only possible for a planet in an orbit with e > 
(Tab.HJ. 



6.5.3. Capture in "negative" resonances 

We now look at "negative" resonances, that is, resonances 
with p < (TabQ}. When s < 90°, we saw that the equilibrium 
rotation rate is always positive (Eqs|92land ll08l ) and negative 
spin states could not be attained. However, for retrograde plan- 
ets (whose obliquities are higher than 90°), the equilibrium ro- 
tation rate is negative (Eq J108t . Since capture in a "positive" 
resonance can hardly occur for these values of the obliquity 
(Eq |109l l, the rotation rate will decrease until the planet starts to 
rotate in the opposite direction. Once in this situation, the evo- 
lution for negative rotation rates can be depicted from the posi- 
tive rotations case. Due to the symmetry in the spin equations, 
the couple {-u, s) behaves identically to the couple {w, n - s). 
The capture probabilities can then be obtained directly from the 
positive case: we just have to take into account that the situa- 
tion corresponding to capture from faster spin rates (slower in 
modulus) will now correspond to the capture from slower spin 
rates and vice versa. Contrary to the case with e = 0°, the 1/2 
resonance will now be the first resonance to be encountered, 
followed by the 1 / 1 . 

Capture in "negative" resonances is then a real possibility 
for planets whose evolution leads the spin to the final obliquity 
s = 180° (section EOt . However, notice that, although p < 0, 
in this case the resonces always correspond to prograde rota- 
tion, since e > 90°. These spin-orbit resonances are thus not 
true negative resonances, and that is why we wrote "negative" 
with quotes. Indeed, an observer of the planet today is unable 
to determine if the 3/2 spin-orbit resonance corresponds to the 
spin state (co/n = 3/2, s = 0°) or (aj/n = -3/2, s = 180°). A 
similar result had already been described for the spin of Venus 
dCorreia and Laskarl 120011) . but for a different kind of equilib- 
rium. 

As a consequence, a true negative resonance will correspond 
to p < and e < 90° (or p > and e > 90°). These resonant 
configurations are not impossible to achieve, but the capture 
probability is very low (Figs|7]and|9]l and we did not register a 
single case during our numerical simulations. 
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6.5.4. Core-mantle friction 

According to expression d64l i the core-mantle friction effect 
on the spin has two different contributions, one permanent if 
the obliquity is different from zero (Eql73l and another arising 
only near spin-orbit resonances CEql60b. This last contribution 
was already present for the case e — 0° and therefore examined 
in detail in section l6".4.2l As before, its contribution to the cap- 
ture probabilities for a non-zero obliquity is given by expression 
( 11091 ) with, for weak coupling, 



2c r K nl ll 



K(l +x 2 )Q(e) 
or, for strong coupling, /3 m = fi x , 

g — and ju = 1 + 



and 



1 +q , 



Ap x n 



3kK{\ + x 2 )Q.{e) 



(HO) 



(HI) 



The contribution of the non-resonant core-mantle friction 
term to the capture probabilities is also easy to compute when 
we use the linear approximation given by expression d73T >. In 
absence of any other dissipative effect, according to expression 
d90l ) it is given by 



cap 



<o. 



(112) 



This probability is always negative, which means that P+ = 
for any resonance. Because the spin can only decrease (Eq|92|i, 
we also have P~ = 0. 

In a more general case, tides are also present and must be 
taken into account. The total dissipative torque can then be 
rewritten as 



D(y) 



where 



1 + x 2 
-K——n(e) 



2xE(e) 
l+x 2 



A 



£ = 2pn 



K f (pn) (l-x 2 ) 



2\ 2 
X X 



K£l(e) 



1 + x z 



>0 



(113) 



(114) 



is the ratio between the magnitudes of the core-mantle friction 
and tidal effects. The total capture probability is now obtained 
straightforwardly from expression d90l > as 



P ± =2 

cap 



Pi P( 



2xE(e) 
1 +x 2 



2Aw 



with 



P( 



: P ~ (- , 
P 



(115) 



(116) 



q — and fx = 1. If we want to take into account the res- 
onant contribution from core-mantle friction, we just need to 
modify g and [i according to expressions (111 01 > and ( 111 11 1. For 
a dominating core-mantle friction » 1) the previous expres- 
sion simplifies to expression (II 121) . On the other hand, when 
( — > (weak core-mantle friction effect) we find expression 
( 1 109b . This is always the case when the obliquity is close to 0°, 
90° or 180°, because x = or x = ±1. 



Since f is always a positive quantity, expression ( 11 151 ) shows 
that the capture probability when de-spinning from faster spins 
is always smaller than it would be without the non-resonant 
core-mantle friction = 0). In figure [9] we plot the global 
effect for different spin-orbit resonances and we observe an im- 
portant reduction in the capture probabilities as the obliquity 
increases. 
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Figure 9: Probability of capture in some spin-orbit resonances when de- 
spinning from faster rotation rates for different obliquities when e = 0.206 
and (B - A)/C m = 1.2 X 10~ 4 . We compute the probability when tidal ef- 
fects are considered alone (solid line) and when core-mantle friction is also 
included, using two different effective viscosities: v = 10~ 6 m 2 s~' (dashed 
line) and v = 10~ 4 m 2 s _1 (dotted line). Although the resonant core-mantle cou- 
pling greatly increases the chances of capture in resonance at zero obliquity 
Ccap = 100% for the 2/1 and 3/2 resonances), for high values of the obliquity 
this probability considerably decreases, allowing the rotation of the planet to 
evolve into lower-order equilibrium configurations. 



We then conclude that, contrary to the resonant contribution, 
the non-resonant core-mantle friction effect is responsible for a 
reduction in the probability of capture for any resonance. Al- 
though the resonant core-mantle coupling greatly increases the 
chances of capture in resonance at zero obliquity (P^ ap = 100% 
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Using the dynamical equations derived in the former sections 
we will now simulate the evolution of Mercury's spin. We start 
our integrations shortly after the formation of the Solar System 
and let the planet evolve until present days. This strategy is 
needed because the present spin configuration of Mercury cor- 
responds to a stable equilibrium and it is impossible to remove 
it from this state by simply reversing the time. 

In order to proceed with our study we need to choose a set 
of plausible coefficients for the dissipative models described in 
sections |3] and |U Although the contribution of these effects to 
the dynamical equations is well understood, some of the geo- 
physical parameters intervening are poorly known. We will 
use a set of parameters called the standard model (&2 = 0.4, 
Q = 50, v = 10~ 6 m 2 s~ 1 ), whose choice is described in detail in 



Correi a and Lask ar (2009). 



7.1. Early evolution 



Before looking at the final stages of the spin evolution, where 
capture in spin-orbit resonances may occur, we will study the 
behavior of the spin in the early stages after planetary forma- 
tion. At this epoch Mercury is supposed to rotate much faster 
than today and any o r ientation of its axis in space is allowed 
( Pones and Tremainei 19931; Kokubo and Ida , 20071) . Indeed, 
the Caloris Basin, a large multi-ringed impact structure is es- 
timated to have been formed by the impact of a 150 km object 
about 3.85 Gyr ago, at the end of the period know n as late heavy 
bombardment dMurdinl bOOOt IStrom et al.ll2008l) . 



Because Mercury is believed to spin rapidly at the begin- 
ning of its evolution {w » «), a constant Q model ( section l3~2l 
seems to be the best choice for tidal evolution. However, for the 
present slow rotation (u> ~ n), a viscous tidal model is the most 
appropriate. In our simulations we then interpolate between the 
two models, adopting Q — 20 for a fast rotating planet with a 
constant dissipation model, and Q — 50 in the limit of slow ro- 
tations with a viscous dissipation model. The transition is done 
at a) ~ 10 n. 




Figure 10: Obliquity evolution with the rotation rate for an initial rotation pe- 
riod of about 32 h (<x), = 65 n). We observe two distinct behaviors: for initial 
obliquities smaller than 175° the obliquity is brought to zero. For higher initial 
obliquities the final obliquity is 180° and the final rotation rate is negative. 



We follow the spin evolution for an initial rotation period of 
about 32 h (a> = 65 n) and different initial obliquities span- 
ning from 0° to 180° (FigfTUI). Two distinct behaviors are ob- 
served: for initial obliquities smaller than about 175° the obliq- 
uity is brought to zero. For higher initial obliquities the final 
obliquity is 180° and the final rotation rate is negative. As 
discussed in section 16.21 there are only two final possibilities 
for the obliquity. The bifurcation in these two distinct evo- 
lutions is provoked by core-mantle friction. In the fast rotat- 
ing regime (to » n), this dissipative effect can be neglected 
(Eql52l and tidal effects drive the obliquity to the equilibrium 
value e - 67° (Eql38l. This is why initial obliquities lower 
than 67° increase while higher decrease. Once in the slow rota- 
tion regime (to ~ «), the tidal equilibrium obliquity will move 
toward zero (Eq|93l. However, in this regime, core-mantle fric- 
tion becomes stronger than tidal effects and drives the final 
obliquity evolution. When core-mantle friction becomes dom- 
inant, if the obliquity is still higher than 90° it will then end at 
180° (Eq|9"Tl). In absence of core-mantle friction the obliquity 
would always end at zero degrees. 

For a stronger core-mantle friction effect, the global picture 
in figure [10] would not change much, only the critical initial 
obliquity that triggers the two distinct final evolutions would 
be lower than 175°. For instance, using v = 10" 4 m 2 s _1 this 
threshold drops to 170°. It is also important to note that the 
critical obliquities increase for faster initial rotation rates and 
decrease for slower ones. 

7.2. Final evolution 

We now look in detail at the behavior of the spin at the final 
stages of its evolution. Here the planet will encounter several 
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spin-orbit resonances where the rotation rate may be trapped. 
Since we are considering the effect of core-mantle friction, we 
need to take into account its resonant contribution to the rota- 
tion rate (Eql60l). as discussed in section l3T2l With the presently 
known geophysical parameters of Mercury, using the core vis- 
cosity v = 10~ 2 m 2 s~', we compute k 2 //3 ~ 10~ 2 <K 1, the k 2 //3 
ratio being smaller for lower viscosity values. Thus, we con- 
clude that we can use the weak coupling approximation (Eqf77T> 
for Mercury, and the evolution of its rotation rate can be de- 
scribed by expression (1102l i. Nevertheless, in our simulations 
we will integrate simultaneously the mantle (Eqf54b and the 
core (EqEoTl rotation rates. 

As we just saw in previous section, the final evolution of the 
obliquity is either along 0° or 180°. We will then first consider 
the simplified situation where the obliquity already achieved 
one of the final states (s = 0°) when the rotation rate enters the 
zone of spin-orbit resonances. This will allow us better to com- 
pare our results with those from previous studies, for which the 
obliquity was always held fixed at zero degrees. Because there 
is no guarantee that at the time of the first resonance crossing 
the obliquity is already close to zero (FigfTOll and since the cap- 
ture probabilities also change with the obliquity (section [675b . 
we will then look at the evolution of the rotation rate when the 
obliquity is still varying. 

7.2.7. Cases = 0° 

Once the obliquity reaches 0° (or 180°), the non-resonant ef- 
fect of core-mantle friction vanishes (Eqs l73l9Tl i. Tidal dissi- 
pation will then drive the rotation rate of the planet towards 
a limit equilibrium value cj e depending on the eccentricity e 
and on the mean motion n (Eq|96b. In a circular orbit (<? = 0) 
this equilibrium coincides with synchronization (a> e /n = 1), 
while the equilibrium rotation rate aj e /n - 1.5 is achieved for 
<?3/2 = 0.284927 (Fig|3]l. For the present value of Mercury's 
eccentricity (e = 0.206) we have Lu> e /n = 1.25685. Thus, 
when Mercury is de-spinning from faster rotation rates it will 
encounter all spin-orbit resonances with p > 3/2, and when 
the spin is increasing from lower values it can be captured in 
resonances with p < 1 / 1 . In absence of planetary perturba- 
tions, the eccentricity remains constant and each resonance is 
crossed only once. In order to estimate numerically the proba- 
bility of capture we kept the initial rotation period of 32 h, 
and ran 2000 simulations for an initial obliquity of 0° with only 
slightly different initial libration phase angles. Since the reso- 
nant part of the core-mantle friction contribution CEqj60b mod- 
ifies the probability of capture (section l6~4l . we performed our 
experiments first in absence of core-mantle friction, and then 
including this effect with v = 10~ 6 m 2 s _I . Results are shown in 
Table |2] 

In absence of core-mantle friction and using the present ec- 
centricity of Mercury (<? = 0.206), the probability of cap- 
ture in the 3/2 spin-orbit res onance was numerica ll y estim ated 
to be 7.2%. In their work, IGoldreich and Peald (119661) esti- 
mated analytically the same probability to be P3/2 = 6.7% 
(Eq llOOl i. With the update d value of the mome nt of inertia 
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Table 2: Capture probabilities in several spin-orbit resonances (in percent). 
The first column (Pj ap ) refers to the theoretical estimation given by expres- 
sion U00K while the next column (num.) refers to the estimation obtained 
running a numerical simulation with 2000 different initial conditions. We used 
(B - A)/C m = 1.2 X 10 -4 , and a constant eccentricity e = 0.206. 



with our numerical simulations. The same is also true for all 
the other resonances (Tabf2|l. 

When we add the effect from core-mantle friction, the nu- 
merical simulations are still in a good agreement with the theo- 
retical estimations given by expression ( llOOI i. showing that they 
can be used to forecast the behavior of numerical simulations. 
In this situation the capture probabilities considerably increase 
for all spin-orbit resonances. In particular, capture probability 
in the 2/1 resonance also becomes 100%, preventing a subse- 
quent evolution to the 3/2 resonance. This behavior was al- 
ready expected from our analysis in section 16.4.21 and it is in 
conf ormity with the resul ts from IGoldreich and Pealel d 19671) 
and IPeale and Bossl ( 1977 ). that is, when the effect from core- 
mantle friction is considered, the probabilities of capture are 
greatly enhanced. 



(B - A)/C m - 1.2 x 10~ 4 (lAnderson et aU 11987). this proba 



bility is increased to 7.7%, which is in a satisfactory agreement 



7.2.2. Case s + 0° 

The initial obliquity of Mercury is unknown, since a small 
number of large impacts at the end of the formation process will 
not average away and may change the obliquity of the planet 
jDones and Tremaine , 1993 ; iKokubo and Idal 12007 ). In addi- 
tion, even for initial low obliquities, during the first stages of 
the evolution, the strong tidal effects acting on Mercury tend to 
increase the obliquity f FigfTOb. Thus, when the planet arrives in 
the slow rotation regime (a> ~ n), for which resonance crossing 
occurs, it is almost certain that the obliquity is higher than zero. 

According to expression (1108b as the obliquity increases, the 
equilibrium rotation rate decreases (Fig|6]l allowing the spin to 
evolve into spin-orbit resonances lower than the synchronous 
rotation (p = 1). This possibility is also enhanced by the fact 
that the capture probability in positive resonances (p > 0) is 
also reduced for high obliquities (Figs|7] and [9). Thus, when 
considering a non-zero obliquity we may expect the planet to 
evolve into resonant configurations with p < 3/2. 

To test these scenarios, we repeated the previous 2000 nu- 
merical experiences using several different initial obliquities 
with v = 10~ 6 m 2 s~'. In Table [3] we report the distribution of 
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the different final spins obtained. As expected, we observe a 
significant modification in the number of captures for all reso- 
nances, and for retrograde planets these captures can occur in 
lower-order resonances than the 3/2. 

For initial obliquities lower than the critical obliquity (so w 
175° for v = 10~ 6 m 2 s _1 ) the capture probability in all reso- 
nances is reduced according to expression (II 15b . Indeed, when 
the resonances are crossed with non-zero obliquity, not only the 
capture probability resulting from tides is smaller, but there is 
also a reducing contribution from the core-mantle friction effect 
(Figj9j». For instance, when the initial obliquity is sq = 90°, the 
2/1 resonance is crossed with an obliquity around s - 50°, and 
capture in the 3/2 resonance becomes possible, while it was not 
for s — 0° (Fig|9]i. The fraction of captures in this resonance 
presents some variations between 1% and 7%, because it de- 
pends on the probability of not being captured in higher-order 
resonances (Tabf3]l. 

For initial obliquities higher than 170° (but still lower than 
the critical obliquity 175°), we observe that capture in reso- 
nances lower than the 3/2 also occurred. The reason is that, as 
core-mantle friction decreases the obliquity, it also decreases 
the rotation rate following expression (l92l . When the obliquity 
reaches zero, the rotation rate may be lower than the equilib- 
rium rotation rate determined by tides (Eq|96l>, i.e., < a> < 
E(e)n, with £(0.206) = 1.25685. Then, tides alone will in- 
crease the spin toward the equilibrium position, and if co is be- 
low the spin-orbit resonances 1/2 or 1/1, capture in these res- 
onances can occur. Notice however, that in this case the cap- 
ture probabilities are given by the expression of P™, since the 
planet crosses the resonance when the spin is increasing from 
lower rotation rates. This is why capture in the 1/2 resonance is 
possible (Fig|5j. 

For initial obliquities higher than the critical value, the obliq- 
uity of the planet is always higher than 90° when the resonances 
are crossed (FigfTOt.Thus. the probability of capture in all "pos- 
itive" resonances (p > 0) is very small (Fig|9]l. In this situa- 
tion the equilibrium rotation rate given by tides will be set at 
oj = -1.25685 n (Eq ll08b . A planet with a> > will con- 
tinue to decelerate, skip all the "positive" resonances, reverse 
its rotation direction, and accelerate its spin rate again. It will 
then sequentially cross resonances 1/2 and 1/1 until it reaches 
a) = -1.25685 n, if not captured before in one of those two 
resonances. The capture probabilities in these two "negative" 
resonances are the same as for the positive resonances 1/2 and 
1/1 for a planet with s — 0° when increasing its spin from 
lower values (section l6.5.31 >. In the case with v = 10 -6 m 2 s , 
we have a 31.9% chance of capture in the 1/2 resonance and 
100% in the 1/1, that is, all the simulations that avoided the 1/2 
resonance were then trapped in the 1/1 (Tab|2j. 

8. Conclusions 

In the present work we derived a formalism to describe the 
complete evolution of the spin of a terrestrial planet like Mer- 
cury under the effect of strong solar tides and core-mantle fric- 
tion. The inclusion of the obliquity in the equations of motion 
allowed us to compute the capture probabilities in resonance 



for any initial spin value. Even though zero obliquity is the fi- 
nal equilibrium resulting from dissipative effects, for many sets 
of initial conditions it is possible that Mercury encounters a res- 
onance with a large obliquity. As we increase the obliquity of 
the planet, the probability of capture in resonance always de- 
creases, allowing the spin to evolve to unexpected configura- 
tions. 

The presence of a liquid core and the associated core-mantle 
friction effect also may lead to a peculiar evolution of the spin. 
Indeed, if the obliquity is higher than 90° at the moment this 
effect becomes dominant over tides, the final obliquity will be 
set at 180°. In this case the rotation rate evolves into negative 
values, and the final spin is always prograde. 

Another important consequence of core-mantle friction is to 
produce significant modifications of the capture probabilities in 
resonance. This effect can be decomposed in two, one resulting 
from the libration of the mantle near spin-orbit resonances (res- 
onant contribution) and another resulting from different orienta- 
tions of the core and the mantle spin vectors when the obliquity 
is not zero (non-resonant contribution). This last effect leads to 
a decrease in the capture probabilities in resonance, while the 
first one increases those chances. 

Finally, we performed some numerical simulations for Mer- 
cury, starting with a fast rotating planet and different obliquity 
values. This allowed us to study the final distribution possibil- 
ities for the rotation rate. Higher-order resonances than the 3/2 
were already expected, since Mercury had to cross them when 
de-spinning from faster rotation rates. However, we could also 
observe lower-order resonant configurations such as the syn- 
chronous or the 1/2 spin-orbit resonance. For retrograde plan- 
ets, "negative" resonances (p < 0) were also observed, but also 
corresponding to prograde final states, since s = 180°. The 
present formalism and results should apply more generally to 
any extrasolar planet or satellite with a core and whose evolu- 
tion led to cross spin-orbit resonances. 

The synchronous resonance can also be achieved for zero 
obliquity if the chaotic evolution of the eccentricity is taken 
into account, when very low values of the eccentricity desta- 
bilize higher-order resonanc es and drive the pl a net's r otatio n 
towards the 1/1 resonance dCorreia and Laskarl l2004 [2009). 
This cannot occur for the present value of the eccentricity 
(<? = 0.206). In a forthcoming study we will include the con- 
tribution of planetary perturbations, which will require mas- 
sive numerical simulations. This will add the contribution 
of the variation of e ccentr icity, already taken into account by 
Correia and Laskar d2009h . but in addition, using the formal- 
ism that has been developed in the present work, we will be 
able also to take into account the large chaotic variations of 
the planet's obliquity resu lting from planetary perturbations 
dLaskar and RobutelLfl993T) . 
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Table 3: Capture probabilities (in percent) in several spin-orbit resonances for different initial obliquities. Tidal effects and full core-mantle friction are included 
with v = 10~ 6 m 2 s~'. We performed a numerical simulation with 2000 close initial conditions with (B - A)/C,„ = 1.2 X 10~ 4 and e = 0.206. We see that crossing a 
resonance with an obliquity different from zero substantially modifies the chances of being captured in a specific spin-orbit resonance. (*): these final equilibria are 
achieved for a final obliquity of 1 80° and a negative rotation rate. 



A. The mean potential energy 



where </> = w + ifr. Let 



In section [2721 we eliminate the fast angles 6 and v from the 
potential energy 1A by expanding the true anomaly v in series of 
the mean anomaly M and then by taking the average of 1A over 9 
and M. However, resonant terms with argument (G—pM) appear 
in the expression of the potential energy (Eqs. [5] |6]) that must 
be taken into account to the mean potential energy 1A (Eq.QjJ. 
Here we give the derivation of the amplitudes f3 x and /?, as well 
as the respective phase angles <p x and </>,- 

Let U r be the resonant part of the potential energy (Eq.|5]) 



(117) 



where /? is given by expression (TT21 . We can rewrite H r as the 
real part of 



%<U r 
~C~ 



/ a 

-M- 



2(1 -x 2 )e l2s 



+ (x+l) 2 e i(2e - 2w) + (x-l) 2 e i(2e+2w) 
and, averaging over 9 and w using expressions ^ and (0, 



(118) 



C 



-13 J] {2(l-x 2 )G(p, e )e i2 



+ H(p,e)[(x+l) 2 e -^e i2(f> - pM) 

+(x-l) 2 e l2 V 2(fl+ " M) ]} {119) 



g = j(l - x 2 )G(p, e) and /i ± = -{x ± l) 2 H(±p, e) . (120) 
4 8 

Retaining only the the resonant terms with argument (9 - pM) 
in expression dl 19K we rewrite expression (II 17b simply as: 



^ = -y8[g + /i + e-^+/i-e i2 *]e 



\2{e-pM) 



(121) 



where 



.0 



= (g + h + + h~) - 4g {h + + h~) sin 2 (f> 



-4h + h- sin 2 : 



and 



tan 2<f> x 



(h + - h ) sin^ 



(122) 



(123) 



g + (h + + h ) cos 2(p 

If x — 1, we have /3 X = /3H(p,e) and </> A - (f>. Notice that if 
g, h* > (which is often the case), we have 



/3 x </3(g + h + +h~), 
and also that the average value over <p is: 



■ g 2 + (h + ) 2 + (h-f 



(124) 



(125) 



The amplitude a r and the phase angle <f> r can be similarly 
obtained from expression ( I12U . Indeed, from expression (151 
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we have: 
cIe 
dt 



and 



1 



8wsine 



d_ d_ 

C d0 + dt// 



c 



L\ gr + h> -^ + h -e^\ 



i sin e e 



\2{6-pM) 



JC126) 



0) 



where 

^ = -^G(p,e) and hf = j(x ± l)H(±p,e) . (127) 

The expressions for at a r and r are respectively given by the 
expressions for/?, and <p x (Eqs. ll22lll23t . where the quantities 
g and /z* are respectively replaced by g r and hf. 



B. Determination of the expressions for the functions £l(e) 
and N(e) 



In section 13.2.21 we wrote the spin equations of motion 
(Eqs. [33]) under the effect of tides for a viscous dissipation 
mod el. For e = 0, these eq uati ons were a lready obtained 
by iGoldreich and Peald d 19661) and lSul dl98ll) . Here we will 
derive those expressions for any v alue of the eccentric ity and 
obliquity. We already used them in lLevrard et al.1 (120071) . with- 
out demonstration. Using the same notation as in sections 12.11 
and l3.11 we rewrite the tidal potential (Eq.l22l as: 



Gm Q R 
r r 



-P 2 (cosS) 



(128) 



with cos S = f f ', the prime ' referring to the interacting body. 
Assuming that both interacting and perturbing body are in the 
same orbital plane, we can write 



cos S 



(\+x) 2 



„2 r 



cos(w - w' - 6 + ff) 



cos(w -w' +6-0) 



cos(vv - w') 



+ cos(w + w') cos(0 - ff) - 1 



(J29) 



where W — m' + ifr + v' is the true longitude of date. It is now 
easy to evaluate the contributions to the spin variations using 
equations (l24l l. For the variation of the rotation rate we obtain: 



dL 
dt 



3Gm'm R 5 d cos S 

■ k.2 5—5 cos S 

80 



rV 3 



(130) 



Let At be the time delay between the perturbation and the 
planet's response. Then, assuming At small and the interact- 
ing body the same as the perturbing one (m' = m Q ), we write 
(Mignard. ^979lll980h : 



ff = 0(t - At) =* 0(t) 



d6^ 
dt 



-At m 6 - coAt 



(131) 



v = v(t - At) =! v(0 - —At = v-n% Vl - e 2 Af . (132) 
dt r 1 

Substituting the above expressions (1131b and (11321 ) into expres- 
sion ( 11301 ) for the rotation rate, we have to first order in At: 



dL 
dt 



3Gm 2 B R 5 



1 + X l— XT 



2 2 

,2 



■ cos(2w) 



- [x^j Vl -e 2 \n 



At. (133) 



Using ajr = ( 1 + e cos v)/( 1 - e 2 ) in the previous expression and 
averaging it over the mean anomaly M and the longitude of the 
perihelion nr, we finally get: 



dL _ Gm%R 5 3k 2 
dt ' a 6 Q 

where Q~ l = nAt, 



(l+x 



2\ 



Q(e) xN(e) 

n 



(134) 



Q(e) = 



1 + 3e 2 + 3e 4 /8 



and 



N(e) = 



l + 15e 2 /2+45e 4 /8+5e 6 /16 



(135) 



(136) 
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dimensionless parameter 
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body rigidity 
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y 


kinematic viscosity 


I4T1 


Ve 


body viscosity 


M 


$ 


internal structure factor 






longitude of the perihelion 
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symbol 


Designation 


Eq. 


P 


Mercury's mean density 


1311 


Q 


dimensionless parameter 


|100| 


cr 


tidal frequency 


1231 


T a 


time constant for damping body tides 


l30l 


Tb 


time constant for damping body tides 


[301 




libration phase (0 = m + i//) 
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libration phase 
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(f>x 


libration phase 


1 1 T5 I 
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t 


general precession angle 


fT4l 




general precession (i]/ x = <p x + \j/ cos e) 


Pel 

I65I 


a> 


mantle's rotation rate 


I43I 


co c 


core's rotation rate 


[391 


k 


projection of a> c over k 


1551 




initial rotation rate 


1921 




equilibrium rotation rate 


|o/;l 




precession angular velocity 


M 


Q(e) 


tidal eccentricity function 
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